Selection of the scaling solution in a cluster coalescence model 
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The scaling properties of the cluster size distribution of a system of diffusing clusters is studied in 
terms of a simple kinetic mean field model, ft is shown that a one parameter family of mathematically 
valid scaling solutions exists. Despite this, the kinetics reaches a unique scaling solution independent 
of initial conditions. This selected scaling solution is marginally physical; i.e., it is the borderline 
solution between the unphysical and physical branches of the family of solutions. 
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Phenomena of diffusion, aggregation and coalescence 
of clusters of particles occur in various scientific areas 
. Of particular interest to us is the relevance of these 
phenomena to the physics of thin films. In recent years, 
material scientists have been able to observe || clusters 
of atoms or vacancies diffusing on thin film surfaces. The 
diffusion and coalescence of these clusters affect the mor- 
phology of the film, which in turn affects its electrical 
properties and the possibilities to use it in the fabrica- 
tion of electronic devices. 

In many systems, the time-dependent cluster-size dis- 
tribution exhibits a scaling behavior, which has been the 
subject of several existing studies pL||. The theoretical 
studies start from a kinetic model of the diffusing clus- 
ters, and then determine the scaling exponents by making 
a scaling ansatz. Finally, the scaling function is found by 
solving the kinetic model numerically starting from some 
initial cluster size distribution. 

The purpose of this work is to investigate the proper- 
ties of the scaling function. In particular, the questions 
answered in this Letter are: Does the scaling function 
depend on initial conditions? Starting from a scaling 
ansatz, is there a unique solution for the scaling func- 
tion? The answers to both questions turn out to be no. 
We show, on the one hand, that simulations of a cluster 
coalescence model starting from different initial cluster 
size distributions approach the same scaling function in 
the long time limit. On the other hand, solving the equa- 
tion for the scaling function that results from a scaling 
ansatz, we obtain a one parameter family of scaling func- 
tions. This implies that out of all the solutions for the 
scaling function, one is selected dynamically. A selection 
criterion is proposed and verified numerically. The appli- 
cability of this selection criterion to other models, which 
exhibit a scaling behavior, will be studied in future work. 

The simplest model that describes diffusion and coa- 
lescence of clusters is the model of Meakin [0 , where D- 
dimensional spherical clusters diffuse in a d-dimensional 
space with diffusion constants, T>(s), that depend on the 
cluster size, s (the number of particles in the cluster). 
When two clusters with s% and S2 particles touch, they 
merge irreversibly and form a new spherical cluster of 
size s = Si + s 2 . Thus, the cluster size distribution de- 



pends on time, and the average cluster size growth is 
limited only by the finite number of particles in the sys- 
tem. Following Meakin, we assume an algebraic decay of 
the cluster diffusion constant with size: 



V(s) 



(1) 



In the present work we use D = d = 2. This is 
the case of two-dimensional clusters diffusing on a two- 
dimensional surface, appropriate for diffusion of islands 
(or voids) of atoms on the surface of a thin film, for ex- 
ample. In this case, it is possible to show that the cluster 
diffusion constant indeed decays algebraicly with cluster 
size. The different values of the exponent £ correspond 
to various microscopic mechanisms responsible for island 
diffusion ||. For example, £ = 3/2 when mass trans- 
port on the surface is dominated by diffusion of atoms 
along island boundaries. If, however, mass transport is 
dominated by exchange of atoms between islands and 
nearby terraces, one obtains £=1 or £ = 1/2. When 
this exchange process is fast compared with the diffusion 
of atoms on the terraces, £ = 1. The case £ = 1/2 is 
obtained in the opposite limit, when surface diffusion is 
fast compared with the exchange kinetics. 

To analyze the scaling behavior of such a system, we 
adopt a mean field approximation due to Smoluchowski 
M. In this approach, the time-dependent density of clus- 
ters of size s, p(s,t), is assumed to obey the following 
rate equation: 
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[D(s)+V( S ')]p(s,t)p(s',t) 



(2) 



The r.h.s. of this equation consists of a gain term due 
to the coalescence of clusters of sizes s' and s — s' into a 
cluster of size s, and a loss term due to the coalescence 
of clusters of size s with clusters of all sizes. 

Simulations of various versions of this model have 
shown (see, e.g. Q) that, in the long time limit, the den- 
sity of clusters obeys the scaling relation 



p(s,t) = t 



' f [w> 



f(t) ~ t 



(3) 
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where s(t) is the time dependent average cluster size, 
and the exponents a and (3 depend on the form of the 
diffusion constant: a — 2/(^+1) and f3 = + 1). 
These values of the exponents can be easily verified by 
substituting the scaling ansatz (||) in Eq. (||) (see below). 

Here we show that the Smoluchowski equation together 
with the scaling ansatz are not sufficient to find the scal- 
ing function /, since they lead to a one parameter family 
of scaling functions. To show this, we start from the 
generalized scaling ansatz 



p(s,t) = g(t)f 



s 



(4) 



We substitute this form of p(s, t) into Eq. (||), and change 
variables from s and t to u = s/s(t) and t. We also re- 
place summation over s by integration over u. This is 
justified, since in the long time limit s(t) is very large. 
Therefore, when s changes by one, the change in u is 
much smaller than one. These algebraic manipulations 
together with the functional form of the diffusion con- 
stant, V(s) — K/ , lead to the following equation: 



g'(t) 



/(«) 



uf(u) 
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(5) 



where g'(t) and f'(u) are the first derivatives of g{t) and 
f(u), respectively. 

Now we take into account the conservation of the total 
number of particles in the clusters, 0: Y] . sp(s, t) = 0, 
where does not depend on time. Using the scaling 
ansatz (0), we obtain a relation between g(t) and s(t): 



<?(*) 







f 



Jo uf{u)du 



(0) 



This result can now be combined with Eq. ([5J) to get the 
following equation for the scaling function /: 



f 



ds 



K@ dt 



u'f(u')du' [2f(u)+uf'(u)} 



—f(u')f(u-u')du' 

fit's. 



(7) 
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f(u)f(u')du' 



The r.h.s. of Eq. (0) depends on it, but not on t. Hence, 
the l.h.s. of the equation cannot depend on time, and we 
can rewrite Eq. (j?]) in the form 



H u'f(u')du'[2f(u) + uf'(u)] 
Jo 



U -Lf(u')f(u-u')du / + 



(8) 
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f{u)f(u')du' , 



where fj, = 4| / KQ is a constant that does not depend 



on t or u. This immediately implies that s ~ i 1 / < ^ +1 - ) , and 
therefore (3 = + 1) as stated above. Since g ~ s~ 2 , 
the value of the exponent a is also confirmed. 
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FIG. 1. Numerical solutions of the equation for the scaling 
function with ( = 1/2. The different curves represent solu- 
tions with various values of fx (the number near the curves): 
(a) n> fio and (b) (i < fj,o- 

In addition, we can now find the scaling function f(u) 
by solving Eq. (|J) with the appropriate boundary condi- 
tion. In this case, it is an integral condition that we 
derive from the definition of the average cluster size, 
s = s p( s j t)/ J27Li p( s i t)- Using the scaling ansatz 



2 



iii this definition, we obtain the integral condition 



f(u)du= / uf{u)du . (9) 
ii Jo 



Before we actually solve for the scaling function /, it 
is important to note that one can multiply any solution 
of Eq. (||) by a constant to get a new solution that satis- 
fies the same integral condition. To eliminate this trivial 
(and physically meaningless) freedom, we enforce another 
integral condition: 



f(u)du = 1 



(10) 



Interestingly, the constant fj, is a free parameter in the 
problem of finding the scaling function as defined above. 
Our aim now is to show that there is a range of values of 
the constant /i, for which there are legitimate solutions 
of Eq. (|J) that obey the integral conditions (||) and (10). 
They are all valid scaling functions for this system. 

We solved the equation for the scaling function numer- 
ically for £ = 1/2, using a method that will be outlined 
elsewhere [pi , and a few examples of solutions are plot- 
ted in Fig. [I]. We found that there is a special value of 
/z, which we denote by /j,q w 1.232, such that for /i > /io 
the solutions are perfectly valid. Three examples of such 
solutions with fi = 1.232, 1.9 and 2.5 are shown in Fig. 
|l|(a). The solutions with /i < /xo, on the other hand, are 
unphysical, since f(u) < in these cases in a range of 
values of u. This can be seen in Fig. 0(b) in the cases 
/i = 0.6, 0.75 and 1.0. In fact, for even smaller values of 
/j,, the scaling function develops oscillations, and its limit- 
ing behavior at large values of u becomes ill-defined. We 
conclude that there is a one parameter family of valid and 
physical scaling functions which correspond to [i > 

Our question now is which of these scaling functions 
are actually reached by the physical system? This is 
of course a question about the kinetics of the system, 
and in order to answer it we should solve the kinetic 
Smoluchowski equation (^) starting from different initial 
cluster size distributions. This turns out to be a diffi- 
cult task that consumes enormous amounts of computer 
time, since the amount of computer time required for ev- 
ery time step is proportional to the square of the number 
of possible values of cluster sizes. 

To circumvent this difficulty, we used a much more ef- 
ficient simulation method. Each simulation started with 
a set of 10 6 -10 7 clusters of various sizes picked according 
to the chosen initial cluster size distribution. At each 
time step a pair of clusters were picked at random. Let 
us denote their sizes by si and S2- The two clusters were 
merged into a single cluster of size S1+S2 with probability 
P{si 1 S2) proportional to T>{s\) + P(sa)- The clock was 
then advanced by 5t ~ 1/N(t) 2 , where N(t) is the total 
number of clusters in the system at time t. This process 
was repeated many times thus leading to an evolution of 



the cluster size distribution. It will be shown elsewhere 
H that the evolution of the average number of clusters of 
size s, induced by this kinetic model, follows the Smolu- 
chowski equation. This procedure is much more efficient 
than a simple integration of Eq. (Q), and the only price 
we have to pay is the statistical error induced by the 
stochastic nature of the model. 
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FIG. 2. (a) A log-log plot of island size distributions re- 
sulting from the simulations with £ = 1/2, at various stages 
of the evolution of the system, (b) The scaled island size dis- 
tributions of (a) form a single curve to a very good accuracy. 

The simulations were carried out for the case £ = 1 /2 
(results for other values of ( will be presented elsewhere 
H ) , starting from several different initial cluster size dis- 
tributions. For example, we started from clusters which 
were all of size 25, and also from a uniform distribution 
of clusters between the sizes 1 and 50. All the simula- 
tions we did exhibited a scaling behavior, and the scaling 
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function was independent of the initial distribution to the 
degree of accuracy of the results. In Fig. ||(a) we show 
the cluster size distributions, at different times, obtained 
from one of the simulations starting from a set of 10 6 
clusters, all of size 1. We see from the figure that the 
most probable cluster size as well as the width of the 
distribution grow with time. Fig. |^(b) shows t a p(s,t) 
as a function of the scaled cluster size u = s/s(t), with 
a = 4/3 as deduced from the scaling analysis for £ = 1/2 
(see above). Clearly, there is excellent data collapse and 
all the distributions fall on top of a single curve, which is 
the scaling function f(u) up to a multiplicative constant. 
Examination of the average cluster size as a function of 
time shows that s(t) ~ t@ with (3 = 2/3 as expected from 
the scaling analysis for ( = 1/2 (see above). 
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FIG. 3. The scaling function obtained from the simula- 
tions is shown as full circles. The solid line is the solution of 
the scaling function equation with fj, = /Lto ~ 1.232. 

Our simulation results indicate that the scaling func- 
tion is not sensitive to the initial cluster size distribu- 
tion. This intriguing result implies that one scaling func- 
tion (and a unique value of the parameter fi) is selected 
dynamically out of the one parameter family of possi- 
ble physical scaling functions. Which one is it? What 
is the selected value of /i? The solution to this puz- 
zle can be inferred from Fig. ||, which shows the scaling 
function obtained from simulations (normalized so that 
Jo f( u )du = 1) as full circles, and the scaling function 
for fi = /l«o ~ 1-232 as a solid line. The resemblance of the 
two functions suggests that the kinetically selected scal- 
ing solution is the one at the border between the physical 
and unphysical branches of the family of solutions, i.e., 
the scaling function that corresponds to /i = (j,q. 

One can argue that this selection principle is reason- 
able, since the initial cluster size distribution is cut off 
at some finite size. Assuming that the tail of the distri- 



bution develops gradually, the system will approach the 
physically legitimate distribution with the shortest tail. 
A close examination of the family of solutions leads to 
the conclusion that the scaling function with ji = fiQ has 
the shortest tail, since it decays to zero faster than all 
the scaling functions with > jiQ. 

A similar selection principle was proposed by Stavans 
et al. and Segel et al. || in relation with the selec- 
tion of the steady state in the coarsening of cells in two- 
dimensional soap froths, and by Maggs et al. jlO| in con- 
nection with steady state length distributions of living 
polymers. This suggests that the selection principle pro- 
posed here may apply to a class of problems. In future 
work, we intend to investigate its applicability in various 
non-equilibrium systems that exhibit a scaling behavior. 

This research was supported by grant No. 95-00268 
from the United States-Israel Binational Science Foun- 
dation (BSF), Jerusalem, Israel. I would like to thank 
E. Domany, T. L. Einstein, Y. Manassen, D. Mukamel, 
Z. Olami, J. Stavans and J. D. Weeks for discussions. 
D. Kandel is the incumbent of the Ruth Epstein Recu 
Career Development Chair. 



[1] E-mail address: fekandel@weizmann.weizmann.ac.il. 

[2] See, for example, B. J. Mason, The Physics of Clouds, 
(Oxsford Univ. Press, Oxsford, 1957); N. A. Fuchs and 
S. G. Sutugin, Highly Dispersed Aerosols, (Ann Arbor 
Science Publications, Ann Arbor, 1970); D. Toussaint 
and F. Wilczek, J. Chem. Phys. 78, 2642 (1983). 

[3] J.-M. Wen, S.-L. Chang, J. W. Burnett, J. W. Evans 
and P. A. Thiel, Phys. Rev. Lett. 73, 2591 (1994); J.-M. 
Wen, J. W. Evans, M. C. Bartelt, J. W. Burnett and P. A. 
Thiel, Phys. Rev. Lett. 76, 652 (1996); K. Morgenstern, 
G. Rosenfeld, B. Poelsema and G. Comsa, Phys. Rev. 
Lett. 74, 2058 (1995); E. Ter-Ovanesyan, Y. Manassen, 
N. Ramesh Rao and Z. Olami, J. Vac. Sci. Tech. B (in 
press) . 

[4] P. Meakin, Physica A 165, 1 (1980). 

[5] D. S. Sholl and R. T. Skodje, Physica A 231, 631 (1996). 

[6] S. V. Khare, N. C. Bartelt and T. L. Einstein, Phys. Rev. 

Lett. 75, 2148 (1995); S. V. Khare and T. L. Einstein, 

Phys. Rev. B 54, 11752 (1996). 
[7] M. v. Smoluchowski, Phys. Z. 17, 585 (1916); M. v. 

Smoluchowski, Z. Phys. Chem. 92, 129 (1917). 
[8] D. Kandel, unpublished. 

[9] J. Stavans, E. Domany and D. Mukamel, Europhys. Lett. 
15, 479 (1991); D. Segel, D. Mukamel, O. Krichevsky and 
J. Stavans, Phys. Rev. E 47, 812 (1993). 
[10] A. C. Maggs, D. Mukamel and C. A. Pillet, Phys. Rev. 
E 50, 774 (1994). 



4 



